Experimentos com pêndulos¶
In [1]:
using Plots
using Images
using ImageDraw
using VideoIO
using Base64
using CSV
using DataFrames
using DifferentialEquations.OrdinaryDiffEq
using Optim
using StatsBase
In [2]:
include("scripts/DisplayMovie.jl")
include("scripts/ImageTracking.jl")
Out[2]:
In [3]:
filename = "img/pendulo_70cm_1_reduzido_curto.mp4"
vd = VideoIO.load(filename)
Out[3]:
 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
(a vector displayed as a row to save space)
In [4]:
typeof(copy(vd[1]))
Out[4]:
In [5]:
println("Duração do vÃdeo: $(VideoIO.get_duration(filename)) s")
println("Número de quadros: $(length(vd))")
println("Número médio de quadros por segundo: $(length(vd)/VideoIO.get_duration(filename))")
In [6]:
# video completo
filename_completo = "img/pendulo_70cm_1_reduzido.mov"
vd_completo = VideoIO.load(filename_completo)
println("""Duração do vÃdeo: $(VideoIO.get_duration(filename_completo)) s""")
println("Número de quadros: $(length(vd_completo))")
println("""Número médio de quadros por segundo: $(length(vd_completo)/VideoIO.get_duration(filename_completo))""")
In [7]:
vd[1]
Out[7]:
In [8]:
DisplayMovie.display_movie("img/pendulo_70cm_1_reduzido_curto.mp4")
In [9]:
background = convert.(RGB{Float16}, vd[1])
Out[9]:
In [10]:
convert.(RGB{Float16}, vd[54])
Out[10]:
In [11]:
convert.(Gray{Float16}, abs.(convert.(RGB{Float16}, vd[54]) - background))
Out[11]:
In [12]:
threshold = 0.1
mask = convert.(
Gray{Float16},
abs.(
convert.(RGB{Float16}, vd[54]) -
background
)
) .> threshold
Out[12]:
In [13]:
convert.(Gray{Float16}, mask)
Out[13]:
In [14]:
labeled_components = label_components(mask)
Out[14]:
In [15]:
@which label_components(mask)
Out[15]:
In [16]:
?label_components
Out[16]:
In [17]:
colorset = distinguishable_colors(maximum(labeled_components)+1, colorant"black", dropseed = false)
Out[17]:
In [18]:
labeled_components
Out[18]:
In [19]:
coloredmask = map( n -> colorset[n+1], labeled_components)
Out[19]:
In [20]:
maximum(labeled_components)
Out[20]:
In [21]:
components_location = ImageTracking.locate_components(labeled_components)
min_area = 1
filter!(p -> p.area ≥ min_area, components_location)
length(components_location)
Out[21]:
In [22]:
components_location = ImageTracking.locate_components(labeled_components)
min_area = 10
filter!(p -> p.area ≥ min_area, components_location)
length(components_location)
Out[22]:
In [23]:
vd54 = copy(vd[54])
for p in components_location
draw!(vd54, Polygon(ImageTracking.RectanglePoints(p)), colorant"red")
end
vd54
Out[23]:
In [24]:
components_location = ImageTracking.locate_components(labeled_components)
min_area = 100
filter!(p -> p.area ≥ min_area, components_location)
Out[24]:
In [25]:
vd54 = copy(vd[54])
for p in components_location
draw!(vd54, Polygon(ImageTracking.RectanglePoints(p)), colorant"red")
end
vd54
Out[25]:
In [26]:
components_location = ImageTracking.locate_components(labeled_components)
min_area = 200
filter!(p -> p.area ≥ min_area, components_location)
Out[26]:
In [27]:
p = components_location[1]
draw(vd[54], Polygon(ImageTracking.RectanglePoints(p)), colorant"red")
Out[27]:
In [28]:
filename = "img/pendulo_70cm_1_reduzido_curto.mp4"
tracks_70cm1 = ImageTracking.find_tracks(filename)
Out[28]:
In [29]:
ImageTracking.tracks_on_video(filename, tracks_70cm1)
In [30]:
DisplayMovie.display_movie("img/pendulo_70cm_1_reduzido_curto_tracked.mp4")
In [31]:
tracks_70cm1
Out[31]:
In [32]:
tracks_70cm1[1].nspan
Out[32]:
In [33]:
tracks_70cm1[1].path
Out[33]:
In [34]:
tracks_70cm1[1].path[1]
Out[34]:
In [35]:
plt = plot(title = "Paths of $(length(tracks_70cm1)) track(s)", titlefont = 10,
xlims = (1,854), ylims = (1, 480),
size = (854, 480), aspect = :equal, legend=:topright)
for (n, tr) in enumerate(tracks_70cm1)
scatter!(plt, getfield.(tr.path, :x), 481 .- getfield.(tr.path, :y), label="track $n")
end
display(plt)
In [36]:
plt = plot(title = "Coordinate x of $(length(tracks_70cm1)) tracks", titlefont = 10, legend=:topleft)
for (n, tr) in enumerate(tracks_70cm1)
scatter!(plt, tr.nspan, getfield.(tr.path, :x), label="track $n")
end
display(plt)
In [37]:
plt = plot(title = "Coordinate y of $(length(tracks_70cm1)) tracks", titlefont = 10, legend=:topright)
for (n, tr) in enumerate(tracks_70cm1)
scatter!(plt, tr.nspan, - getfield.(tr.path, :y), label="track $n")
end
display(plt)
In [38]:
size(vd[1])
Out[38]:
In [39]:
data_x = getfield.(tracks_70cm1[1].path, :x)
data_y = size(vd[1], 1) + 1 .- getfield.(tracks_70cm1[1].path, :y)
nothing
In [40]:
n_minima = Int[]
for j in 40:length(data_x)-1
if data_x[j-1] > data_x[j] < data_x[j+1]
push!(n_minima, j + first(tracks_70cm1[1].nspan) - 1)
end
end
n_minima
Out[40]:
In [41]:
plt = plot(title = "Coordinate x of $(length(tracks_70cm1)) tracks", titlefont = 10, legend=:topleft)
for (n, tr) in enumerate(tracks_70cm1)
scatter!(plt, tr.nspan, getfield.(tr.path, :x), label="track $n")
end
vline!(plt, n_minima, label = "minima")
display(plt)
In [42]:
(n_minima[2:end] - n_minima[1:end-1]) / 29.5
Out[42]:
In [43]:
g = 9.8 # m/s
l = 0.7 # m
T_p = 2π * √(l / g)
Out[43]:
In [44]:
function medias(a::Real, b::Real)
(a ≥ 0 && b ≥ 0) || throw(ArgumentError("arguments must be nonnegative"))
ma = (a + b)/2
mg = sqrt(a * b)
return ma, mg
end
function agm(a::Real, b::Real; tol::Real = 1e-10, maxiter::Int = 100)
tol > 0 || throw(ArgumentError("tolerance must be positive"))
maxiter > 0 || throw(ArgumentError("maximum number of iterations must be positive"))
y1, y2 = a, b
n = 0
while abs(y1 - y2) > tol && n < maxiter
y1, y2 = medias(y1, y2)
n += 1
end
return (y1 + y2)/2, abs(y1 - y2), n
end
Out[44]:
In [45]:
θ₀ = π / 4
T_p / agm(1, cos(θ₀ / 2))[1]
Out[45]:
Escala¶
- Vamos fazer ajustar um circunferência ao movimento do pêndulo, para descobrir a escala espacial.
In [108]:
error(u) = sum(
abs2,
(data_x[50:end] .- u[1]).^2 + (data_y[50:end] .- u[2]).^2 .- u[3]^2
) / length(data_x)
Out[108]:
In [109]:
center_point = Optim.optimize(error, [400.0, 800.0, 500])
Out[109]:
In [110]:
Optim.minimum(center_point)
Out[110]:
In [49]:
center_nx, center_ny, radius_n = Optim.minimizer(center_point)
Out[49]:
In [105]:
plt = plot(title = "Paths of $(length(tracks_70cm1)) track(s)", titlefont = 10,
xlims = (1,854), ylims = (1, 600),
size = (854, 600), aspect = :equal, legend=:topright)
for (n, tr) in enumerate(tracks_70cm1)
scatter!(plt, data_x, data_y, label="track $n")
end
scatter!(plt, (center_nx, center_ny), markersize = 10)
plot!(plt, 1:854, nx -> center_ny - √(radius_n^2 - (nx - center_nx)^2), linestyle = :dash)
display(plt)
In [51]:
scale = radius_n / l # points per meter
Out[51]:
In [52]:
data_scaled_x = (data_x .- center_nx) / scale
data_scaled_y = (data_y .- center_ny) / scale
Out[52]:
In [53]:
function dudt_pendulum!(du, u, p, t)
θ, ω = u
l, g, α = p
du[1] = ω
du[2] = - (g / l) * sin(θ) - α * ω
end
Out[53]:
In [54]:
g = 9.8 # 9,8 m/s^2
l = 0.7 # 70 cm
α = 0
p = [l, g, α]
θ₀ = π/6
ω₀ = 0.0
u₀ = [θ₀, ω₀]
tspan = (0.0, 5.0) # até 5 segundos
prob_pendulum = ODEProblem(dudt_pendulum!, uâ‚€, tspan, p)
Out[54]:
In [55]:
sol = solve(prob_pendulum, Tsit5())
Out[55]:
In [ ]:
In [56]:
display(plot(sol, label = ["θ(t)" "ω(t)"]))
display(plot(sol, vars = 1, label = "θ(t)"))
In [57]:
n0 = 34
plt = plot(sol, vars = 1, label = "θ(t)")
scatter!(plt, (0:length(data_x)-n0) ./ 29.5, data_scaled_x[n0:end], label="tracked")
display(plt)
In [58]:
df = CSV.read("data/pendulo_70cm_reduzido.csv", DataFrame)
Out[58]:
In [115]:
plt = plot(title = "Coordinate x", titlefont = 10, legend=:topleft)
scatter!(plt, df.nf, df.nx, label = nothing)
display(plt)
plt = plot(title = "Coordinate y", titlefont = 10, legend=:topleft)
scatter!(plt, df.nf, size(vd[1], 1) .+ 1 .- df.ny, label = nothing)
display(plt)
In [114]:
data2_scaled_t = (df.nf[n0:end] .- df.nf[n0]) / 29.97
data2_scaled_x = (df.nx[n0:end] .- center_nx) / scale
data2_scaled_y = (size(vd[1], 1) .+ 1 .- df.ny[n0:end] .- center_ny ) / scale
display(scatter(data2_scaled_t, data2_scaled_x))
display(scatter(data2_scaled_t, data2_scaled_y))
In [62]:
data_x2 = getfield.(tracks2[1].path, :x)
n_minima2 = Int[]
for j in 40:length(data_x2)-1
if data_x2[j-1] > data_x2[j] < data_x2[j+1]
push!(n_minima2, j + first(tracks2[1].nspan) - 1)
end
end
n_minima2
In [116]:
n_minima2 = Int[]
for j in 40:length(df.nx)-1
if df.nx[j-1] > df.nx[j] < df.nx[j+1]
push!(n_minima2, j + df.nf[1] - 1)
end
end
n_minima2
Out[116]:
In [117]:
(n_minima2[2:end] - n_minima2[1:end-1])/30
Out[117]:
In [118]:
sum(n_minima2[2:end] - n_minima2[1:end-1])/7/30
Out[118]:
In [122]:
plt = plot(title = "Coordinate x", titlefont = 10, legend=:topleft)
scatter!(plt, df.nf, df.nx, label = "data")
vline!(plt, n_minima2, label = "minima")
display(plt)
In [66]:
g = 9.8 # 9,8 m/s^2
l = 0.7 # 70 cm
α = 0
p = [l, g, α]
θ₀ = π/6
ω₀ = 0.0
u₀ = [θ₀, ω₀]
tspan2 = (0.0, 16.45)
prob_pendulum2 = ODEProblem(dudt_pendulum!, uâ‚€, tspan2, p)
Out[66]:
In [67]:
VideoIO.get_duration("img/pendulo_70cm_1_reduzido.mov")
Out[67]:
In [68]:
θ₀ = π / 6
α = 0.0 #-0.1
remake(prob_pendulum2; u0 = [θ₀, 0.0], p = [l, g, α])
sol2 = solve(prob_pendulum2, Tsit5(); saveat = 1/29.87)
data2_scaled_x[1:10] .- sol2(data2_scaled_t)[1,1:10]
#length(sol[1,:])
#sol(data2_scaled_t)[1,:]
sol(data2_scaled_t)
Out[68]:
In [69]:
plot(data2_scaled_t, sol(data2_scaled_t)[1,:], label = "ODE", ylims=(-1,1))
scatter!(data2_scaled_t, data2_scaled_x, label = "data")
Out[69]:
In [70]:
function ddu_pendulum!(ddu, du, u, p, t)
θ = u[1]
ω = du[1]
l, g, α = p
ddu[1] = - (g / l) * sin(θ) - α * ω
end
Out[70]:
In [71]:
θ₀ = π / 6
ω₀ = 0.0
α = 0.05
l = 0.74
p = [l, g, α]
prob_pendulum_2nd = SecondOrderODEProblem(ddu_pendulum!, [ω₀], [θ₀], tspan2, p)
sol_2nd = solve(prob_pendulum_2nd)
sol_2nd.retcode
Out[71]:
In [72]:
plot(data2_scaled_t, sol_2nd(data2_scaled_t)[2,:], label = "ODE", ylims=(-1,1))
scatter!(data2_scaled_t, data2_scaled_x, label = "data")
Out[72]:
In [137]:
prob_rmk = remake(prob_pendulum_2nd; u0 = [θ₀], p = [l, g, α])
sol = solve(prob_rmk)
In [94]:
function objective(β)
g = 9.8
ω₀ = 0.0
θ₀, l, α = β
p = [l, g, α]
tspan2 = (0.0, 16.45)
prob = remake(prob_pendulum_2nd, u0 = [θ₀, 0.0], p = [l, g, α])
#sol = solve(prob)
prob = SecondOrderODEProblem(ddu_pendulum!, [ω₀], [θ₀], tspan2, p)
sol = solve(prob)
return mse(data2_scaled_x, sol(data2_scaled_t)[2,:])
end
Out[94]:
In [95]:
res = Optim.optimize(objective, [Ï€ / 3, 0.7, 0.0])
Out[95]:
In [83]:
θ₀, l, α = Optim.minimizer(res)
Out[83]:
In [76]:
π / θ₀
Out[76]:
In [77]:
θ₀ - π / 5.6529
Out[77]:
In [78]:
remake(prob_pendulum_2nd, u0 = [θ₀, 0.0], p = [l, g, α])
Out[78]:
In [79]:
sol_2nd = solve(prob_pendulum_2nd)
sol_2nd.retcode
Out[79]:
In [80]:
plot(data2_scaled_t, sol_2nd(data2_scaled_t)[2,:], label = "ODE", ylims=(-1,1))
scatter!(data2_scaled_t, data2_scaled_x, label = "data")
Out[80]: